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TECHNICAL MEMORANDUM 


THE NASA/MSFC GLOBAL REFERENCE ATMOSPHERIC MODEL - 

1990 VERSION (GRAM-90) 

Part I: Technical/Users Manual 


1. INTRODUCTION 


1.1 Background 

In response to needs for engineering reference model atmospheres of wide scope and appli- 
cation, Georgia Tech developed, under NASA sponsorship, a Global Reference Atmospheric Model 
(GRAM) with latitude, longitude, and monthly variations over a height range of 0 to 700 km 
(Justus et al., 1974a,b, 1975, 1976). The original GRAM version has undergone a series of 
improvements over the years (Justus et al., 1980, 1986, 1988). This report describes additions and 
improvements to this model whereby extensive new data sets from the middle atmosphere (20 to 
120 km) have been incorporated. This new version is called GRAM-90. This section describes the 
basic GRAM-90 version and section 2 provides a more detailed description of the GRAM-90 
program. Section 3 presents some sample results for the middle atmosphere height region. Sections 
4 and 5 are user’s and programmer’s manuals for the GRAM-90 program. A complete listing of 
the GRAM-90 code and middle atmosphere input data set is available in part II of this report (con- 
tact Dale Johnson, ES44, NASA Marshall Space Flight Center, Huntsville, AL 35812). 


1.2 Description of the Basic Model 

The 1990 version of GRAM-90 is an amalgamation of two previously existing empirical 
atmospheric models for the low (<25 km) and high (>90 km) atmosphere, with the new empirical 
latitude-longitude dependent model for the middle atmosphere (20 to 120 km). Smooth transition 
between these height regions is insured by fairing techniques in the overlap height regions. In the 
original GRAM version, the high atmospheric region above 1 15 km was simulated by the Jacchia 
(1970, 1971) model. In GRAM-90, a version of the Jacchia model called the Marshall Engineering 
Thermosphere (MET) model is used (Hickey, 1988a,b). The atmospheric region between 25 km 
and 120 km is simulated by a new latitude-longitude dependent empirical model, based on exten- 
sive new data sets which have come primarily from the Middle Atmosphere Program (MAP). 

These data sources and the resultant new data set used in GRAM-90 are described more fully later 
in this report. Between 90 and 120 km, a smooth transition between the middle atmosphere data 
and the MET values is accomplished by a fairing technique. Below 25 km the atmospheric 
parameters are computed by a four-dimensional (4-D) world-wide atmospheric model developed for 
NASA by Allied Research Associates (Fowler and Willard, 1972; Spiegler and Fowler, 1972). 
Between 25 and 30 km, an interpolation scheme is used between the 4-D results and the middle 
atmosphere values. Figure 1 . 1 presents a schematic summary of the atmospheric regions of the 
GRAM-90 program and how they are modeled. 


Height, 


1000 


MET (JACCHIA) MODEL 
(Hickey, 1988a, b) 


120 


FAIRING BETWEEN THE MIDDLE ATMOSPHERE 
DATA AND THE MET VALUES 


I 90 


NEW DATA FROM THE MIDDLE ATMOSPHERE PROGRAM (MAP) 
(References in text) 


30 


INTERPOLATION BETWEEN THE MIDDLE ATMOSPHERE 
DATA AND 4-D VALUES 


25 


SURFACE 


4-D MODEL (Spiegler and Fowler, 1972) 


Figure 1.1. Schematic summary of the atmospheric regions in the GRAM-90 program, and the 
simulation methods used for mean monthly values in each region. 
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1.3 New Zonal Mean Data 


For the middle atmosphere section, described more fully in section 2, monthly zonal mean 
data from six different sources have been merged: 

(1) Groves (1985) zonal mean thermodynamic and zonal wind data for 18 to 80 km. 

(2) Barnett and Comey (1985a) zonal mean thermodynamic and zonal wind data for 20 to 
80 km. 

(3) Fleming et al. (1988) zonal mean thermodynamic and zonal wind data for 20 to 
120 km. 

(4) Zonal mean thermodynamic zonal wind and meridional wind data from Dartt et al. 

(1988) for 20 to 85 km. 

(5) Koshelkov (1985) zonal mean thermodynamic and zonal wind data for Southern Hemi- 
sphere meteorological rocket data (20 to 65 km). 

(6) Zonal mean thermodynamic, zonal wind, and meridional wind data processed from the 
latest (1969-1986) “SUMS” data for the Meteorological Rocket Network (MRN) for the 
Northern Hemisphere (20 to 65 km). 

These data have been processed at 10°-latitude intervals from 80° S. to 80° N. and at 5-km 
height intervals from 20 km upward (to the top of the respective data sets). Values for the poles 
(90° S. and 90° N.) are added by an extrapolation procedure. 

1.4 New Stationary Perturbation Data 

Monthly latitude-longitude cross section data (for revision of the stationary perturbation data) 
have been prepared from three different sources: 

(1) Groves (1985) planetary wave I and 2 amplitudes and phases for thermodynamic data 
over 18- to 80-km height. 

(2) Barnett and Comey (1985b) planetary wave 1 and 2 amplitudes and phases for thermo- 
dynamic data over 20- to 80-km height. 

(3) Latitude-longitude cross sections of thermodynamic, zonal wind, and meridional winds 
from Dartt et al. (1988) data for 20 to 85 km. 

These data have been reduced to relative perturbations about the zonal-mean value (sta- 
tionary perturbations) for the 10° latitude increments from 80° S. to 80° N., 20° longitude intervals 
(180°, 160° W., 140° W., ..., 160° E.), and 5-km height intervals from 20 km upward (to 80 or 
85 km). Stationary perturbations above 85 km were estimated by a vertical extrapolation technique. 
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1.5 New Middle Atmosphere Wind Model 

Observational data for zonal-mean zonal wind (at 10° latitude increments) and latitude- 
longitude variations of zonal and meridional wind (at 10°-latitude by 20°-longitude increment) have 
been added to the data base for the middle atmosphere. These replace the spherical harmonic wind 
model of the previous versions of GRAM. 

Below 25 km and above 90 km, the winds are computed from horizontal pressure gradients, 
estimated by finite differences, from the geostrophic wind relations. Near the equator, an interpola- 
tion scheme is used instead of the geostrophic relation (which approaches infinite values as the lati- 
tude approaches zero). In the height range 25 to 90 km, the newly developed data base of middle 
atmosphere winds is used, at all latitudes. Mean vertical winds, of the order of 1 cm/s or less, are 
evaluated from the slopes of isentropic surfaces and the horizontal advective winds. Wind shear in 
the monthly mean horizontal wind is estimated from horizontal temperature gradients in the regions 
below 25 km and above 90 km. In the 25- to 90-km height range, mean wind shear is computed 
from finite differences of the winds in the new middle atmosphere data base. 


1.6 The Quasi-Biennial Perturbations 

In addition to the monthly mean values of pressure, density, and temperature, two types of 
perturbations are evaluated: quasi-biennial oscillations (QBO) and random perturbations. The QBO 
variations in pressure, density, temperature, and winds, empirically determined to be represented by 
an 870-day-period sinusoidal variation, have amplitudes and phases which vary with height and 
latitude. The QBO amplitudes are primarily significant at low altitudes (—20 to 40 km) for equa- 
torial latitudes, and at higher altitudes (50 to 60 km) for higher latitudes. 


1.7 The Random Perturbations 

For realistic simulation of actual atmospheric parameter values, as they would likely be at 
any given time, random perturbations are also computed and applied as perturbations to the month- 
ly mean values. The random perturbations are evaluated by a simulation technique which uses 
empirical values of magnitudes and scales for the perturbations, in order to generate random pertur- 
bations which have realistic space and time correlations. Complete descriptions of the random 
perturbation model have been presented by Justus et al. (1986, 1988) and Justus (1988). 

For the GRAM-90 random perturbation model, standard deviations about the monthly means 
have been obtained and processed for three data sets: 

(1) Zonal-mean values of standard deviations of thermodynamics, zonal wind, and meridi- 
onal winds from the Dartt et al. (1988) data for 20 to 85 km. 

(2) Zonal-mean values of standard deviations for thermodynamics, zonal wind, and meridi- 
onal wind from the MRN “SUMS” tape (20 to 65 km). 

(3) Monthly standard deviations at 44° N., 6° E. from Chanin et al. (1985) lidar data. 
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These data have been processed as zonal-means of the standard deviations at IO°-Iatitude 
intervals for 80° S. to 80° N. and height intervals of 5 km for 20 km upward. 

For GRAM-90 a new vertical component for the random perturbation winds has also been 
added, based on the turbulence model of Justus et al. (1990). 

1.8 Overview of This Report 

The following sections give a technical description of the GRAM-90, with emphasis on the 
new additions, and new user’s manual descriptions of the program operation aspects of the revised 
model. Section 2 provides a more detailed technical description of the GRAM-90 program. Section 
3 presents some sample results for the new middle atmosphere section and provides comparisons 
with results from a 3-D global circulation model. Section 4 is a user’s manual which provides 
details on how to run GRAM-90. Section 5 is a programmer’s manual which provides more details 
for those wishing to make their own GRAM program adaptations. 


2. TECHNICAL DESCRIPTION OF THE MODEL 
2.1 The Jacchia Section (Above 90 km) 

The Jacchia (1970, 1971) model for the thermosphere and exosphere was originally 
implemented to compute atmospheric density at satellite altitudes. The Jacchia model accounts lor 
temperature and density variations due to solar and geomagnetic activity, diurnal and semiannual 
variations, and seasonal and latitudinal variations. The Jacchia model assumes a uniformly mixed 
composition from sea level to 105 km, with diffusive equilibrium among the constituents (nitrogen, 
oxygen, argon, helium, and hydrogen) above 105 km. Fixed boundary values for temperature and 
density are assumed at 90 km. Alterations, described in Justus et al. (1974a), were made to allow 
atmospheric pressure to be computed from the density and temperature. Geostrophic winds are 
evaluated in the Jacchia section by computing horizontal pressure gradients with successive 
evaluations of the Jacchia model at different latitudes and longitudes. For GRAM-90, the NASA 
MET model (Hickey, 1988a,b) has been implemented as the version of the Jacchia model, for 
altitudes above 120 km. Between 90 and 120 km, a fairing routine, described below, is used to 
insure a smooth transition between the middle atmosphere and the MET model. 


2.2 The 4-D Section (Below 25 km) 

The 4-D atmospheric model, developed by Allied Research Associates (Spiegler and Fowler, 
1972) was designed to extract from data bases and interpolate on latitude and longitude, mean 
monthly and daily variance profiles of pressure, density, and temperature, at 1-km intervals from 
the surface to a height of 25 km for any location on the globe. The data bases contain empirically 
determined, atmospheric-parameter profiles at a large array of locations. The Northern Hemisphere 
grid array is equivalent to the NMC grid network. Grids spaced at 5° intervals of latitude and 
longitude are used in the equatorial and Southern Hemisphere regions. 
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Technical changes previously made in the 4-D program were: a modified latitude-longitude 
interpolation method, previously described in Justus et al. (1974a); an adjustment routine to modify 
the variance to comply with the constraints imposed by the perfect gas law (Buell, 1970, 1972); 
and a check routine to determine vertical and horizontal consistency of the 4-D data. 

For the GRAM-90 program, a routine has been incorporated to read the 4-D data by ran- 
dom access from disk files, rather than by serially reading through the 4-D data tape files. The 
random access feature computes the record number in the 4-D File, based on the latitude and 
longitude required. The specific record is then read directly by accessing that record number from 
the data base. This feature should provide a substantial increase in speed of execution for simula- 
tions along trajectories which traverse long distances at altitudes below 30 km (e.g., for the 
National Aero-Space Plane (NASP) and the High Speed Civil Transport (HSCT) programs). 

The method of application of the 4-D model in the GRAM-90 program is as follows: the 
first time that atmospheric parameters at a location below 30 km are required, a set of atmospheric 
profiles of monthly mean and daily variances of pressure, density, and temperature are generated at 
a 16-point grid of locations spaced at 5° latitude and longitude intervals (a slightly different grid is 
used near the poles and equator). This grid of profiles, covering 15° x 15° of latitude-longitude, is 
then stored in the computer, and all further atmospheric parameter values in the 0- to 25-km range 
are found by interpolation between locations within this grid. If the trajectory goes outside this grid 
while the height remains below 25 km, the program executes an additional call on the routine 
which sets up the 4-D data grid. 

The location of the grid points to be evaluated is determined dynamically, based on the 
position and direction of travel along the trajectory, when the 4-D grid is first required, by a 
procedure described in Justus et al. (1974a). The 4-D data bases normally contain data for the 
surface to 25 km in 1-km steps. At locations where the surface is at more than 1 km above sea 
level, the surface value will be followed by one or more zero records, and the first nonzero record 
above the surface value will be at the lowest integer km higher than the surface. For example, if 
the surface is at 700 m, then there will be data at surface, I km, 2 km, etc.; but if the surface is 
at 1.3 km, the data will contain the surface, one zero record, 2 km, 3 km, etc. An interpolation 
routine (based on the hydrostatic relation and constant lapse rate altitude segments) is used to fill in 
data between sea level and the first nonzero data above the surface. Interpolation is also used to 
fill in any missing data immediately below the 25-km height. The basic interpolation equations 
were described in Justus et al. (1974a). 


2.3 The Middle Atmosphere Section 

For incorporation into the revised GRAM, the six sets of zonal-mean data have been 
merged into one set covering the 20- to 120-km altitude range. The latitude-longitude cross-section 
(stationary-perturbation) data were also merged into one set covering the 20- to 85-km height 
range. Both of these sets include explicit Southern Hemisphere data for the first time (as opposed 
to using only 6-month displaced and north-south reflected Northern Hemisphere data, as was done 
in the original GRAM version. These data also include, for the first time, zonal-mean and latitude- 
longitude cross-section data for zonal and meridional wind components. These replace the previous 
spherical harmonic wind model. 
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The starting point for the middle atmosphere (20- to 120-km) section of GRAM-90 is the 
zonal-mean data set, derived for the data sources discussed in section 1.3. Tabulations of the 
middle atmosphere data base are at intervals of 5 km in height, 10° in latitude ( — 90° to +90°), 
and 1 month in time. Since the original data bases included only the latitude range —80° to +80°, 
polar values were added by the extrapolation relation 


} ! ±90 — y ± go 3' ±70^ • 


(2.1) 


This assumes a parabolic variation between the ±70° and ±80° values across the polar latitude 
±90°. 


Longitudinal variation about the zonal-mean data is provided by revised stationary perturba- 
tions, based on new data sets described in section 1 .4. The stationary perturbation, s y , in parameter 
y is given by 


s y =(y— <y >)/<>’> 


(2.2) 


where <y> represents the zonal-mean value of y (i.e., averaged around a circle of fixed latitude). 
Note that the definition of s y makes it be identically zero at the pole. The stationary perturbation s y 
for parameter y is added to the zonal-mean value <y> = z y to produce the longitudinally variable 
value, y(lat., long.), according to the relation 


v(lat., long.) = z^l+ty) . (2.3) 

The latitude-longitude dependent values computed by equation (2.3) are used as the monthly mean 
values for the altitude range 30 to 90 km and as the basis for interpolation or fairing over 25 to 
30 km and 90 to 120 km. 


2.4 Interpolation and Fairing 

The 4-D data are available in the data bases at 1-km height intervals and at 5°x5° latitude- 
longitude grids in the southern and equatorial areas and at the NMC grid locations in the Northern 
Hemisphere. NMC grid profiles are always converted (by interpolation) to 5°x5° grids before 
interpolation to the trajectory locations. The general interpolation requirements for the 4-D section 
are height interpolation over 1 km and latitude-longitude interpolation over a 5°x5° square grid. 

The middle atmosphere zonal-mean data are tabulated at 5-km height intervals and 10°- 
latitude intervals. Interpolation is required between these tabulated locations. The stationary per- 
turbations are evaluated at 10°-latitude and 20°-longitude intervals and at 5-km height intervals 
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between 30- and 90-km altitudes. Interpolation between these tabulated locations is also required. 
For values between 25 and 30 km, interpolation between the 4-D data and middle atmosphere data 
is required. The interpolations are always carried out in the program by doing the latitude (zonal- 
mean) or latitude-longitude (4-D or stationary perturbations) interpolation first, and then doing the 
height interpolation. 

The MET model can be evaluated at any height above 90 km and at any latitude and 
longitude, so no interpolation is required. However, between 90 and 120 km, there is overlap 
between the middle atmosphere data and the MET model, so a fairing procedure is used to effect a 
smooth transition between the middle atmosphere data at 90 km and the MET values at 120 km. 

The method used to interpolate pressure, density, and temperature over a height interval 
between heights z, and z 2 is to assume a linear variation of the temperature and of the logarithm of 
the density. The latitude interpolation for the middle atmosphere zonal-mean data is done by 
assuming linear variation between the latitudes <j>i and <j> 2 (which are at A<()= 10° apart). 
Two-dimensional latitude-longitude interpolation between a square or rectangular array of positions 
at latitudes c})| and <j >2 and west longitudes \| and \ 2 is done by the relation 


F(cJ>,\) = Fo + (F l -F 0 )8<f, + (F 2 -F 0 )8\ + (F 3 -F l F 2 + F 0 )8<j> bk , (2.4) 


where 84> is (cj> — 4> i )/(<f >2 — 4*0 and 5\ is (A. — \|)/(\ 2 — \|). 

To accomplish smooth transition between the middle atmosphere values at 90 km and the 
MET values at 120 km, a fairing technique is used. This fairing technique was described in Justus 
et al. (1974a). The fairing is done only at the altitudes 95, 100, 105, 110, and 1 15 km, i.e., 
heights for which there are middle atmosphere values in the data base. Linear interpolation is then 
used to fill in intervening heights, as discussed in the height interpolation section above. 

Interpolation of the random perturbation magnitudes is done linearly on the variance (cr 2 ) 
rather than linearly on the magnitude (a). This is because the perfect gas law adjustment equations 
(Buell, 1970, 1972) are nearly linear in the variances. Thus, once variances have been adjusted for 
the perfect gas law constraint, their adjustment would tend to be preserved by linear interpolation 
on variances, not magnitudes. 


2.5 Geostrophic Winds 

The eastward (i.e., blowing toward the east) wind component u and northward component v 
can be evaluated from the geostrophic wind equations 


= —(1/p/) dp/dy 

(2.5) 

= (1/p f) dp/d.x , 

(2.6) 
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where p is the density, / is the Coriolis parameter (2 fi sin <j>) and dp/dx and dpi by are the east- 
ward and northward components of the horizontal pressure gradient. For evaluation in the model, 
the pressure gradient terms must be approximated by finite differences. Geostrophic wind values 
are computed in the 4-D height range, by finite differencing of the 4-D pressure data. 

Since density appears in the denominator of equations (2.5) and (2.6), the geostrophic wind 
components u g and v g can become unreasonably large at high altitudes (low densities). Physically 
this condition is prevented from occurring because of the effects of molecular viscosity, which 
become important at high altitudes. Therefore, winds in the upper (MET) height range are 
computed by viscous-modified geostrophic relations 


/ (j u g k v g ) 

(2.7) 

f- +k z 


f(fv g + t u g ) 

(2.8) 

f 2 + k 2 



where k — \x/(pL 2 ), p. is the coefficient of viscosity, and L is a viscous damping scale 
(d 2 u/dz 2 « u/L 2 ). When viscosity is small {k « f), equations (2.7) and (2.8) become equivalent to 
the geostrophic relations (u - u g \ v - v g ). At high altitudes, where viscosity effects are large 
(k » f), equations (2.7) and (2.8) become 


u — ~ (L 2 /\l)bpldx 

(2.9) 

V s* -(L 2 /p.)bp/dy . 

(2.10) 


Thus, as the effects of viscosity go from small to large, the winds computed by equations (2.7) 
and (2.8) go from being parallel to the isobars to becoming parallel to the pressure gradient. 


2.6 Thermal Wind Shear 

The wind shear components dUg/dz and dvg/dz in the 0- to 25-km height range are evaluated 
by the thermal wind equations 


du^dz = - (g/JT) bTIdy 

(2.11) 

dVgldz = ( g/JT) dT/dx , 

(2.12) 
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which is the usual form, leaving off a correction term in dT/dz, which is normally small. The hori- 
zontal temperature gradient terms are estimated by finite differences in a similar manner to the 
pressure gradient components in equations (2.5) and (2.6). 

Thermal wind shears are also computed in the MET height range by relations derived from 
differentiation of equations (2.7) and (2.8) 


du/dz 


f(f dug/dz-k dvgfdz) 

f + k 2 


(2.13) 


f(f dVg/dz + k dUg/dz ) 

dXl< ' L f 1 + k 2 


(2.14) 


with the geostrophic wind shears coming from the thermal wind relations (equations (2.11) and 

( 2 . 12 )). 

Since the ordinary geostrophic winds are inversely proportional to the coriolis parameter/ 
(which goes to zero at the equator), these relations give unrealistically large winds at low latitudes. 
To overcome this problem, interpolation between about + 15° and — 15° latitude is used in 
GRAM-90. This interpolation limit, the “minimum geostrophic latitude,” is specified in the 
program input and can be selected as any value between 5° and 18°. 


2.7 Mean Vertical Winds 

GRAM-90 also evaluates mean vertical winds from the slope of isentropic surfaces. On such 
surfaces, the entropy function i|j is 


4* = C p T + gz + (w 2 + v , 2 )/2 = const. 


(2.15) 


Therefore, on isentropic surfaces 


dty/dt + w<h|i/<)x + v<h|//dy + U’di^/dz = 0 , 


(2.16) 


and, if di \tldt is assumed zero, the vertical wind vv can be solved for as 


VV = — [wdvJ//djr + vdl|;/dy]/(dl|//dz) . 


(2.17) 
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By differentiation of equation (2.15), with the assumption that u and v are the geostrophic winds u g 
and v g , and that du/dz and dv/dz are given by the thermal wind relations, equation (2.17) becomes 


w = -C p [u g (dT/dx) + v g (dT/dy)]/{g + C p (dmz) 

+ (g/jT) [v g (dT/d.x) — u g (dT/dy)]} . 


(2.18) 


Mean vertical winds evaluated by equation (2.18) are generally less than 1 cm/s, and hence are 
realistic values for the large scale mean vertical winds affecting mean meridional circulation. 


2.8 The Quasi-Biennial Perturbations 

In the original GRAM, MRN data from 1964 to 1969 were used to evaluate quasi-biennial 
amplitudes and phases in the height range 25 to 65 km. The quasi-biennial period which produced 
minimum variance, when simultaneously evaluating the annual, semiannual, and quasi-biennial van 
ation, was found to be 870 days. For a later GRAM version, the harmonic analysis was done the 
same’ way with MRN data for 1970 to 1972 added to the original data base. Again, the 870-day 
period was found to produce minimum variance for the QBO winds, while a 900-day period did 
slightly better for the thermodynamic variables. In order to retain a single period, the original 
870-day period was chosen as still the preferable value overall. This model for the QBO perturba- 
tions is retained in the GRAM-90 version as well. 


2.9 The Random Perturbation Model 

GRAM-90 uses a simple, first-order, autoregressive model to compute a perturbation value 
at each new trajectory position from the perturbation value at the previous position. The effects of 
variation of the mean values and standard deviations about the mean are accounted for in the calcu- 
lation of the perturbations. The process of computing the perturbations is exemplified by the 
computation of a new density perturbation pi at the new position (coordinates x 2 , y 2 , and z 2 ), 
where the mean value is p? and the standard deviation about the mean is a 2 . The density perturba- 
tion pi at the previous position (coordinates X], y i, and Z|) has already been computed (or given as 
a starting seed value if position 1 is the first position in the trajectory). The mean density value at 
position 1 is p,, and the standard deviation is a,. The new perturbation value p 2 is then computed 

from 


p 2 = Ap 2 p(/f>i + B q\ , 


where A and B are coefficients which force the proper correlation R between the successive per 
turbations pf and pi, and <7, is a random number (mean = 0, standard deviation - 1) selecte 
from a Gaussian (normal) distribution. 
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If the magnitude of the separation between the positions is r, where 


r = [(.v 2 -.v,) 2 + (.V 2 -yi) 2 + (r 2 -r|) 2 J l - ° , (2.20) 

then the correlation R is assumed to be 

R = <p 2 pi>/(o- 2 (j,) = e~' L , (2.21) 

where the angle brackets denote averaging, and L is the correlation integral scale length. The 
coefficients A and B in equation (2.19) are given by 


A = R a^/cr | 

B = o- 2 [] -/?-]' 2 . 


( 2 . 22 ) 


In the actual perturbation computations, two forms of equations (2.19) and (2.21) are used, one for 
small scale (turbulence, gravity wave, etc.) perturbations and one for large scale (tides, planetary 
waves, etc.) perturbations. 

Calculations tor temperature, pressure, and wind perturbations are done by more compli- 
cated expressions, which include additional terms required to preserve density-temperature cross 
correlations and density-wind cross correlations. More complete descriptions of the random per- 
turbation model, including recent modifications, are given by Justus et al. (1980, 1986, 1988) and 
Justus (1988). 


3. SAMPLE RESULTS 

3.1 Example GRAM-90 Results and Comparison With GRAM-88 

Figures 3.1 and 3.2 show height-latitude cross sections of zonally averaged temperature for 
the GRAM-90 version for the months of January and July. The new zonal-mean temperature con- 
tours are similar to the 1988 GRAM version (GRAM-88) results, except near 70° N. to 90° N. (or 
70 S. to 90 S.) and 70- to 75-km altitude. This difference is largely due to unrealistic values of 
temperature extrapolations to the poles in the earlier GRAM version (see figures 10 and 1 1 of 
Justus and Roper, 1987). 

Figures 3.3 and 3.4 similarly compare zonal-mean density in GRAM-90 for January and 
July (expressed as percent deviation from the U.S. 1976 Standard Atmosphere density). The largest 
differences between GRAM-88 and GRAM-90 density values are near 80 km in the polar regions, 
with January values having larger negative deviations near the North Pole for GRAM-88 and larger 
positive deviations near the South Pole for GRAM-90. 


12 



Figures 3.5 and 3.6 show the zonal-mean zonal wind for GRAM-90 data for January and 

July. 


Magnitudes of the total standard deviations for the month of January are shown for GRAM- 
90 temperature in figure 3.7, for density in figure 3.8, for zonal wind in figure 3.9, and for 
meridional wind in figure 3.10. The standard deviations for temperature and density are quite 
similar for GRAM-88 and GRAM-90. The standard deviations for zonal wind are somewhat larger 
for GRAM-88 near 50° N. and 50-km altitude than for GRAM-90. 

The longitudinal variability for temperature in GRAM-90 is illustrated by figures 3.1 1 and 
3.12, which show latitude-longitude cross sections of temperature at 30 km (fig. 3.1 1) and 80 km 
(fig. 3.12) altitudes, with temperature expressed as percentage deviation from the zonal-mean 
temperature at the respective heights and latitudes. 

Comparable GRAM-90 data on longitudinal variability of density are shown in figures 3.13 
and 3.14 for the month of January, also expressed as percent deviations from the zonal-mean 
value. 


Longitudinal variations of zonal and meridional wind in GRAM-90 are shown in figures 
3.15 to 3.18. Values are in meters per second. 


3.2 Comparison of GRAM-90 Results With 3-D SCM Output 

As noted in section 1 , the data for GRAM-90 was compiled from a variety of sources, 
including both rocket and satellite data in the 20- to 120-km height range. These data sources 
represent averages, for each month, over a period of several years, the longest period being the 
MRN data (1969 to 1986). Four of the data sets (Groves, 1985; Barnett and Comey, 1985a,b; 
Fleming et ai., 1988; and Dartt et al., 1988) are all derived from the same basic set of upper 
atmospheric satellite observations. The excellent consistency among these data sets, and with the 
combined averages for the GRAM-90, are illustrated by comparison of figures 3.1 to 3.18 with the 
extensive plots provided by Groves (1985) and Fleming et al. (1988). 

In this section, a comparison between the empirically-based GRAM-90 and a physically- 
based. 3-D stratospheric circulation model (SCM) is presented. For purposes of comparison with 
the GRAM-90 results, a 3-D stratospheric dynamical/chemical circulation model was run for an 
annual simulation of the global atmosphere to a height of —85 km above the surface. The model 
used was an enlarged version of the quasi-geostrophic. spectral model of Cunnold et al. (1975) 
which includes the relevant ozone chemistry in interaction with the dynamical and thermodynamical 
components. A detailed description of the enhanced dynamical/chemical model is contained in 
Alyea et al. (1985). Briefly, however, the model was formulated with 32 vertical levels spaced in 
constant increments of the logarithm of pressure (the top being along the 0.00348-mb pressure level 
or at a height of —85.4 km) and was truncated in its horizontal spectral resolution after planetary 
wave number 18 (using a triangular truncation scheme). The principal forecast equations represent 
the quasi-geostrophic vorticity and the ozone mixing ratio. Temperature is obtained from a form of 
the "thermal wind" relationship, and vertical motion is computed at each time step from the 
diagnostic 3-D "w-equation." 
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At the Earth's surface, topographic data in spectral form is used to determine the vertical 
motion from the horizontal wind field. Tropospheric heat forcing is specified using a simple 
symmetric (about the Equator) model based on climatological observations (see Cunnold et al., 
1975). At the stratospheric and mesospheric levels, solar heating is determined by the ozone dis- 
tribution, and the position of the Sun while ultraviolet (UV) heating/cooling is approximated using 
a linearization involving an “equilibrium” temperature field toward which the UV heating/cooling 
term attempts to adjust (see Cunnold et al., 1975). 

For the comparison runs, the model was initialized with output from a previous seven-wave 
version of the model (Cunnold et al., 1975) and integrated forward in time for 1 and one-half 
years. These initial conditions resulted from a more than 3-year run of the earlier model. For 
comparison with the GRAM-90 results, monthly mean temperature and zonal wind fields from the 
last July and January runs of the SCM will be used. 

Figures 3.19 and 3.20 show the SCM results for the July zonal-mean temperature (K) and 
zonal-mean zonal wind (m/s) fields, respectively. These are to be compared with figure 3.2 for the 
GRAM-90 July zonal-mean zonal temperature cross section and figure 3.6 for the GRAM-90 July 
zonal-mean zonal winds. As seen from these figures, there is a great similarity between the SCM 
and GRAM-90 representations (but note that the SCM figures top off at 80-km altitude). In the 
wind fields, the polar night jets in the two models are represented by westerly wind maximums 
with speeds in excess of 90 m/s (greater than 1 10 m/s for the SCM) centered on ~50- to 60-km 
altitude and —40° S. to 50° S. latitude. The summer-hemisphere, easterly, upper-level maximums 
are somewhat higher in altitude than the westerly jets, located at a height of approximately 70 km 
in the —40° N. to 50° N. latitude belt but, as expected, contain maximum wind speeds that are 
lower than those of the polar night westerlies. This is particularly true for the SCM results in 
which the easterly winds attain maximum speeds that are less than 40 m/s. In both cases, the wind 
fields undergo an approximately 150-m/s change in passing from the polar night westerly jet to the 
summer hemisphere easterly maximum, but the SCM maintains an approximately 20-m/s stronger 
westerly component throughout. This feature is linked, through thermal wind considerations, to the 
somewhat stronger horizontal temperature gradients of the winter (southern) hemisphere found in 
the SCM below about 50 to 60 km than for the observationally based GRAM-90. 

During January the two models exhibit greater differences. Figure 3.21 (temperature) and 
figure 3.22 (zonal wind) represent the SCM zonal-mean cross sections for January and are to be 
compared respectively to figures 3.1 and 3.5, representing the GRAM-90 results. In these cases not 
only are the westerly and easterly maximum wind speeds different, but their locations also differ 
significantly. For example, the GRAM-90 wind maximums are both located at about 65 km above 
the surface and within the 30° to 40° latitude belts of their respective hemispheres. The two 
maximum speeds are similar, in the 50- to 60-m/s range. The SCM cross sections are not nearly so 
antisymmetrical about the Equator. Again, these differences can be seen in the temperature cross 
sections as well as in the wind data. In particular, the horizontal temperature gradient in the lower 
stratosphere (—20- to 35-km altitude) of the Northern Hemisphere SCM cross section (fig. 3.21) 
shows a fairly large warm patch between —35° N. to 50° N. latitude. This feature is nearly 
undetectable in the GRAM-90 results and no doubt accounts, through thermal wind considerations, 
for the more northerly (60° N.) and lower (50-km) altitude of the SCM January polar night jet than 
in the GRAM-90. Other aspects of the January SCM and GRAM-90 comparisons show less 
deviation. 
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Part of the reason for the observed differences in the SCM results and the GRAM-90 data is 
due to the effects of interannual variability. For example, winters at higher latitudes are affected by 
varying intensities of stratospheric warming conditions. Near-equatorial latitudes are affected by 
quasi-biennial oscillations. For the GRAM-90 data base, the effects of interannual variability are 
averaged out because of the use of the multiyear data bases. The SCM results shown here are for a 
single year of model run. 
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Figure 3.2. GRAM-90 height-latitude cross section for July temperature (K). 
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Figure 3.3. GRAM-90 height-latitude cross section for January density 
(percent from U.S. 1976 standard). 




Figure 3.4. GRAM-90 height-latitude cross section for July density 
(percent from U.S. 1976 standard). 
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1-90 latitude-longitude cross section of January density deviations 
from zonal-mean (percent) at 30-km altitude. 
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Figure 3.14. GRAM-90 latitude-longitude cross-section of January density 
from zonal-mean (percent) at 80-km altitude. 
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from zonal-mean (m/s) at 80-km altitude. 
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4. USER’S MANUAL 


4.1 General Program Overview 

The GRAM-90 version of the program is designed to produce atmospheric parameter values 
either along a linear path (to be called a profile), with automatically stepped constant height, lati- 
tude, and longitude increments, or along any set of connected positions (to be called a trajectory) 
which must be input individually into the program. 

There are three general types of input to the GRAM program: (I) a set of three lines, called 
the initial data, which contain the values of the program options, the initial position, the profile 
increments, and other information required before the calculations are begun; (2) a data file 
(SCIDAT9) containing parameter values for the middle atmosphere zonal-mean data, the stationary 
perturbations (deviations from the zonal-means, to produce longitude varying monthly means), and 
random and quasi-biennial perturbation parameter values; and (3) the 4-D data files with one data 
file for each month, containing profiles of monthly mean pressure, density, temperature, and their 
variances from the surface to 25 km, for the entire globe. If it is desired to compute atmospheric 
parameters along a trajectory instead of a linear profile, then a fourth type of data — the trajectory 
times and positions — must be input. 

In terms of program function, the major elements of the GRAM program are the main seg- 
ment (GRAM), the subroutine SCIMOD, which is a driver for all of the atmospheric evaluation 
subroutines, and SETUP, a subroutine used to read the SCIDAT9 data file, and load the necessary 
starting conditions for execution. 

Output of the GRAM-90 program consists of monthly mean pressure, density, temperature, 
wind and wind shear, total (mean plus perturbation) values of pressure, density, temperature, and 
winds, as well as perturbation values and magnitudes. 

A complete discussion of the input, output, and program operation characteristics for the 
GRAM-90 program is given in the following sections of the user’s manual. Section 5 is a 
programmer’s manual giving additional program details which might be useful to users wishing to 
modify the program for specific applications. 


4.2 The 4-D Data Files (0 to 25 km) 

The world-wide meteorological data set developed for the 4-D model by Allied Research 
Associates (Spiegler and Fowler, 1972) was originally binary tapes labelled WW1A to WW3A. 
These binary tape files, one per month, have now been converted to ASCII files which can be read 
by random access (or “direct access") I/O operations. 

Within each file are 3,490 records representing the values at individual grid points of lati- 
tude and longitude. These points are grouped into three grid areas: 288 points on the Northern 
Hemisphere Equatorial (EQN) grid; 1 ,977 points on the Northern Hemisphere (National Meteoro- 
logical Center (NMC)j grid; and 1,225 points on the Southern Hemisphere (SH) grid. On the NMC 
grid, the data were computed at NMC points and are stored in the order given by the NMC grid 
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table given in the SCIDAT9 data file. For the other two grid areas, the data are given at 5° lati- 
tude-longitude intersections westward from the Greenwich Meridian to 5° E. The EQN grid covers 
the latitudes from 0° to 15° N. with points occurring in the following order: 1-^4 = long. 0, lat. 0, 
5. 10, 15; 5-8 = long. 5° W., lat. 0, 5, 10, 15; ...285-288 = long. 5° E., lat. 0, 5, 10, 15. The 
SH grid contains all data from 5° S. to the South Pole as follows; I = South Pole, 2-18 = long. 
0, lat. -5 to -85; 19-35 = long. 5° W.. lat. -5 to -85; ...1,209-1,225 = long. 5° E., lat. -5 to 
-85. It should be noted that the South Pole is given only once, as the first point of the SH dataset. 

Each record consists of 1,268 ASCII characters (formatted 208F6.2, 2F4.2, and 3F4.0). The 
first 208 parameter values contain the thermodynamic data for a latitude-longitude grid point. The 
last five variable values are identifiers. The data are arranged by altitude level for each parameter; 
thus, the first 26 values contain the pressure means from the surface to 25 km, and the next 26 
values contain the pressure variances for the same levels. This pattern continues for the 26 levels 
of temperature means and variances, moisture means and variances, and density means and vari- 
ances. Temperatures are in K, pressures are in mb (multiplied by 100 for use in GRAM-90 to give 
N/m 2 ), and densities are in g/m 3 (divided by 1 ,000 for use in GRAM-90 to give kg/m 3 ). GRAM-90 
does not use the moisture means and variances. 

Variables 209 and 210 contain the latitude and longitude of the data grid point (in tenths of 
a degree). The latitude is always positive (since the Southern Hemisphere is identified by grid), 
and the longitude is always west. 

The last three values (variables 211 to 213) contain three identifier parameters. The first of 
these is the homogeneous moisture region in which the point lies, the second is the point number, 
and the third is the month. It should be noted that the 4-D points are numbered within the grid that 
contains them, and not by their location in the file. Thus, the point numbers run from 1 to 288 
(EQN) grid, 1 to 1,977 (NMC grid), and I to 1,225 (SH grid), not from 1 to 3,490. 


4.3 The SCIDAT9 Data File 

This section describes in detail the data contained in the SCIDAT9 data file. 

NMC Grid Data. This data set gives the 4-D Northern Hemisphere point number and the 
dual index for the corresponding NMC location. The NMC grid locations form an octagonal array, 
centered on the North Pole. The points are at square grid locations on the polar projection used for 
the NMC grid. A conversion between the latitude and longitude (treated as polar coordinates on the 
flat NMC grid plane) and the NMC grid indices (treated as Cartesian coordinates on the projection 
plane) is accomplished by a polar-to-Cartesian coordinate transformation, via equations programmed 
into the 4-D model. The NMC grid data in the SCIDAT9 File merely establishes the equivalence 
between the sequential 4-D NMC point number and the 2-D x-y NMC grid point location. The 
NMC grid data constitute the first group of data in the SCIDAT9 file. The NMC grid data file 
contains 396 FORTRAN readable records with code “N” and 15 integers (A2, 1517 format) in each 
record. 


Zonal-Mean Data. The zonal-mean data consists of 12 monthly sets of zonal-mean values 
for pressure, density, temperature, and zonal wind, tabulated at 10° latitude intervals (from -90° to 
+ 90°) and 5-km height increments (from 20 to 120 km) for each month. Prefix codes, ZP, ZD, 
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ZT and ZU indicate pressure, density, temperature, and zonal wind, respectively. Each record c - 
tains the code, the month, the height in km, and the -90°, -80°, .... 80 , 90 latitude values o 
parameter expressed as a four-digit integer, with an exponent common to all of the values in the 
rerappearing at the end of the record. Thus, a value of 2,761 with an exponent at the end of the 
record of -6 would be the same as 2,761 x lOr 6 = 2.761 x 10 \ Pressure data are in units of 
N/m 2 density values are in kg/m 3 , temperatures are in K, and zonal winds are in m/s. The zona 
mean data set contains 1,008 FORTRAN readable records, with the code and 22 integer values in 

each record (format A2, 14, 15, 1916, 14). 

Stationary Perturbations. The stationary perturbations are latitude-longitude dependent 
relative perturbations to be applied to the zonal-mean values. Data for each of FI 2 i months are given 
for both Northern and Southern Hemisphere latitudes. Prefix codes SP, SD, ST. SU and SV md 
cate stationary perturbation values for pressure, density, temperature, zonal (eastward), or mend.- 
onal (northward) wind components, respectively. 

Each record contains the code, (he month, the height in km. the latitude (-80 to +80) in 
degrees, and then 18 values of stationary perturbations, in per mil <«/10> for >he ^emodynanttc 

variables and 0 I m/s for the winds a. longitude 180», 160" W.. 140" W 140 E„ 160 E. 

The monthly mean value y, for parameter y at any latitude and longitude can be competed rom 
the zonal-mean value z,. at the latitude, and the stationary perturbation s, <m per mil) 
tude and longitude, by the relation 


- (I + s y l 1,000) 


(4.1) 


Note that the stationary perturbation values at 90° latitude are always zero. The stationary perturb 
ation data consists of 15,300 FORTRAN readable records, with a code and 21 integer values 

(format A2, 2115). 

The Random Perturbation Data. Random perturbation magnitudes (standard deviations) 
are latitude dependent only. Prefix codes RP, RD, RT, RU, and RV indicate random perturbation 

magnitudes in pressure, density, temperature, zonal wind, and meridional wind .^^foli owed* 
tively. Each random perturbation record has the code, the month, and the height.n km,/o lowed 
by 19 values of random perturbation magnitude, at 10° latitude increments -90 to +90 followed 
hv a common exponent value. These data give the relative standard deviations cr p Ip , a P /p, and 
a y i”en.)Tr use in ,he random perturbation model. The code RU and RV da.a are similar, 
except that the wind perturbations are absolute deviations in m/s, and cover the height range 0 
200 km, whereas the RP, RD, and RT data cover 20 to 200 knv Random perturbation magm u es 
for 0 to 25 altitudes are provided by the 4-D data base for the thermodynamic vanables The 
dom perturbation data consist of 1,596 FORTRAN readable records with a code and 22 integer 
values (A2, 14, 15, 1916, 14) in each record. 

Larae-Scale Fraction Data. From daily difference analysis described in section 2 of 
Justus et al (1980), the fraction of the total variance (a 2 from the random perturbation data) con- 
tained in the large scale perturbations has been determined as a fraction of height and latitude. The 
SCIDAT9 file contains the annual average fraction (expressed as per mil) of total variance con- 
lamed in Ihe large scale. Large-scale and small-scale magnitudes <x,, and a, are computed from the 
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fractional data f L in per mil (code P for pressure, density, and temperature or code PW for winds) 
by the relations 


= (V^/l,000)a r 

(4.2) 

(V 1 — //„ /1 ,000)o- r , 

(4.3) 


where a T is the total perturbation magnitude. The code P and code PW data sets each contain 25 
FORTRAN readable records, with code word P or PW, followed by 17 integer values (A2, 1717) 
on each record. 

Density- Velocity Correlations. Daily difference analysis (Justus et al., 1980) was also 
used to evaluate the cross correlations R up and R vp for use in the velocity perturbation model, 
described in section 2. Both large-scale and small-scale values of the density-velocity correlations 
were evaluated and are given in the SC1DAT9 data file (codes CL and CS) in per mil (i.e., divide 
by 1,000 to get correlations in the range — 1 to +1). 

The code CS and CL data consist of 50 FORTRAN readable records, with code word CS or 
CL, followed by 12 integer values (A2, 1217) in each record. 

The Quasi-Biennial Oscillation (QBO) Data. The QBO data consist of height- and 
latitude-dependent amplitudes and phases for quasi-biennial variations in pressure (QP), density 
(QD), temperature (QT). and eastward and northward wind components (QU and QV, respec- 
tively). The amplitudes of the QBO wind components are in decimeters per second (0. 1 m/s). The 
phases of all of the QBO parameters are measured in days after January 0, 1966, for the 
occurrence of the first maximum value. Since the period of the QBO variations is taken to be 870 
days, the phases can vary from 0 to 870. 

Each QBO data record contains the code, the height in km, the amplitude and phase for 10° 
latitude, the amplitude and phase for 30° latitude, etc., out to the amplitude and phase for 90 lati- 
tude. There are 80 FORTRAN readable records in the QBO data set. Each record contains 1 1 inte- 
ger values, following a code word (QP, QD, QT, QU, or QV) in format A2, 1517. 


4.4 The Initial Input Data 

The initial input data consists of two free-field (no set format with commas after each 
number) lines containing initial position data, program options, and other information required to 
begin computation, plus an optional third free-field line to give initial random perturbation data if 
random perturbations are to be computed, plus an optional set of trajectory position data (followed 
by a backup card), if trajectory positions are to be read in rather than a linear proFile generated 
automatically in the program. Appendix B gives a brief summary of the input characteristics, a 
summary of the data file setup, and some sample input and output for the program. The following 
gives a more detailed description of each program input line. 
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Input Line Number 1 . The first input line, read in by the main program segment GRAM, 
in free-field format, contains the following information. Designation R indicates real quantities, I 
denotes integer quantities. 

1. Initial Height (R), HI: The initial height in km for the beginning point of the profile or 
trajectory. This can be any nonnegative real number. Atmospheric parameters are never evaluated 

at the first position, which is used only to establish the initial conditions. 

2. Initial Latitude (R), PHIl: The latitude of the initial position, in degrees, with southern 
latitudes negative. If the initial latitude, or any subsequent latitude is greater than 90° in absolute 
magnitude, then a transformation is made: 


lat. = ( 1 80° — |lat.|)(lat./|lat.|) 


(4.4) 


long. = long. + 180° . 


(4.5) 


3. Initial West Longitude (R), THET1: The west longitude of the initial position in 
degrees. East longitude can be put in as negative or converted to 0° to 360° west longitude. If 
negative (east) longitudes are input, they are converted to the 0° to 360° scale before being used by 
the program. At any time during the run, if a longitude gets outside the 0° to 360° it is put back 
into that range by adding or subtracting 360°, as necessary. 

4. I0.7-cm Flux (R), F10: The solar 10.7-cm radio noise flux in units of 10~ 22 W /m 2 (the 
normal units for this parameter) at the time for which the atmospheric values are to be computed. 
This factor is used only in the Jacchia (MET) section, so a value of zero can be used on input if 
the height never goes above 90 km. A value of 230 for both design steady state conditions and for 
maximum conditions may be used, or consult MSFC for monthly predictions. 

5. Mean 10.7-cm Flux (R), F10B: The 81 -day mean solar 10.7-cm radio flux. This 
parameter is used in the Jacchia (MET) section to compute the nighttime minimum global exo- 
spheric temperature (equation 14 in Jacchia, 1970). Use zero if the height does not above 90 km. 

A value of 230 may be used for both design steady-state or maximum conditions, or consult MSFC 
for monthly predictions. 

6. Geomagnetic Index (R), AP: The geomagnetic index a p , used to compute a 
geomagnetic correlation to the exospheric temperature, in equation (22) of Jacchia (1970). Use zero 
if the height does not go above 90 km. A design steady-state value of 20.3 and a maximum- 
condition value of 400 may be used for a p , or consult MSFC for monthly predictions. 

7-9. Date (I), MN, IDA, IYR: The date, for the starting time of the trajectory or profile 
evaluation in month/day/two-digit year form, as three integer input values. The day of the month 
and the year have no direct effect on the program calculations, except for the Jacchia section alti- 
tudes or the quasi-biennial oscillation terms. The month is used to establish which middle 
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atmosphere data, stationary perturbation data, and random perturbation data (including large scale 
fractions and velocity-density correlations) to load from the SCIDAT9 data file into the working 
arrays. The program will work more efficiently if multiple trajectories or profiles are evaluated dur- 
ing one run operation and the months are the same. (This avoids repeated look-up of the middle 
atmosphere, stationary perturbation, and random data from the SCIDAT9 file). 

10-12. Greenwich Time (I), IHRO, MINO, ISECO: The Greenwich mean time for the 
starting position in hours, minutes, and seconds as three integer values. Only the Jacchia (MET) 
section is directly affected by the time of day, so unless the height goes above 90 km, the starting 
time would serve merely as a reference parameter for the particular run being done. Greenwich 
time corresponding to a local time of 0900 hours should be used for design steady-state Jacchia 
(MET) section conditions, and for maximum conditions the local time should be taken as 1400 
hours. 


13. Latitude Increment (R), DPHI: If a linear profile is to be generated automatically, 
this is the latitude increment (in degrees) between successive profile positions. The new latitude 
would be the old latitude plus the latitude increment. For a profile with decreasing latitude (going 
southward) the increment must be negative. Use zero if separate trajectory position input is to be 
read in. If a vertical profile (i.e., changing only height) is to be evaluated, then use zero latitude 
increment. 

14. West Longitude Increment (R), DTHET : If a linear profile is to be generated 
automatically, this is the west longitude increment (in degrees) between successive profile posi- 
tions. The new longitude will be the old longitude plus the longitude increment. For a profile 
progressing eastward, use a negative increment. Use zero if separate trajectory position input is to 
be read in. If a vertical profile is to be evaluated, then use zero increment. 

15. Height Increment (R), DH: The height decrease in km between successive positions, 
for an automatically generated linear profile. The profiles are normally generated downward 
(descending height). (New height = old height - the height increment). If an upward generated 
profile is desired, the height increment should be negative. Downward generated profiles will be 
evaluated until the height is incremented to a negative value or until the maximum number ol posi- 
tions (item 16, first line) is exceeded. 

16. Maximum Number of Positions (I), NMAX: The maximum number of profile posi- 
tions to be generated automatically. This does not include the initial position, for which no atmos- 
pheric parameters are evaluated. Use zero if trajectory positions are to be read in. 

17. Time Increment (I), INCT: The time displacement (seconds) between successive 
automatically generated profile positions. This would normally be set to zero, but could be used as 
a counter to be printed out in the time position with the output. For trajectories, the time for each 
position is read in with the position data (see trajectory input section). The hours, minutes, and 
seconds parameters (read in as items 10—12, first line) are updated according to the new time 
generated by the time increment. However, only the elapsed time in seconds is printed out on the 
output. 
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18. Trajectory Option (I), IOPT: This option tells the program whether a trajectory or a 
linear profile is to be evaluated. A value of zero means that a linear profile is to be generated 
automatically from the parameters read on the first line. A value greater than zero means that 
trajectory position data must be read in to determine the positions at which atmospheric parameters 
are to be evaluated. The unit from which the trajectory data are to be read is specified by the 
(nonzero) trajectory option. Thus, if trajectory data are to be read in from unit 5, use a trajectory 
option of 5. 

19. Output Option (I), IOPP: This option tells the program whether or not to produce 
nonprint output of the atmospheric parameters (see the output description section). Nonprint (i.e., 
disk) output is convenient to use as input to plotter programs. A value of zero means no nonprint 
output. A value greater than zero means to output the data on the unit number equal to the output 
option value. 

20. Minimum Geostrophic Latitude (R), GLAT: Below this latitude (in absolute magni- 
tude) winds are computed by latitude-interpolation between the geostrophic values at latitudes 

+ GLAT and —GLAT. Above this latitude, the usual geostrophic relations are used. Values of 
GLAT between 5° and 18° may be input. Values below 5° are set to 5°; values above 18° are set 
to 18°. 


With normal numbers of decimal places and no unnecessary blank spaces, the above 20 
items should fit into one line. However, if they occupy more than the number of columns allowed 
on one line (usually 80) they may be spread out onto two lines if the appropriate rules of free-field 
input are observed. Consult the operations or FORTRAN manual for your particular system for 
characteristics of free-field input. 

Input Line Number 2. The second input line is read in by the subroutine SETUP and 
contains various unit numbers to be used and various options controlling the random and quasi- 
biennial calculations and other functions of the program. The unit numbers are the parameters used 
in read statements in the FORTRAN program to control which file is being read from. The unit 
numbers are required in the input in order to give maximum flexibility in choice of I/O devices for 
the program. All input items on line number 2 are integers. 

1. SCIDAT9 Input Unit, IUS: This is the unit number of the SCIDAT9 file. If the 
SCIDAT9 file has been assigned by control statements, which establish unit 3 for the SCIDAT9 
file, then the SCIDAT9 input unit number should be three on this input line. 

2. 4-D Input Unit, IU4: This is the unit number for the 4-D data file. Any available unit 
number can be used. If the 4-D file containing the January data, DAT4D0I, has been assigned by 
control statements establishing this as unit 4, then the 4-D input unit number is four on line 2. 4-D 
data files must be named DAT4D0I for January, DAT4D02 for February, ..., DAT4DI2 for 
December. 

3. Random Option, IOPR: This option tells the program whether or not to compute ran- 
dom perturbations. If the value is one, random perturbations are computed. If the value is two, 
then random perturbations are not computed. If any values other than one or two are input, the run 
is terminated with a message “ERROR IN SETUP INPUT” and a dump of the parameters most 
recently read in. 
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4. QBO Option, IOPQ: This option tells the program whether or not to compute QBO 
perturbations. If the value is one, QBO perturbations are computed. For two no QBO perturbations 
are computed, and for any other values the ‘ERROR IN SETUP INPUT” and dump of most recent 
parameters read in is given. 

5. First Random Number, NR1: This number is required as a starting parameter for the 
random number generating subroutine RANDOM. Any positive integer 1 to 30,000 can be used. 

Use a value of one for a standard design application run. Provided all other input is the same, a 
given value for the starting random number will always produce the same random perturbation 
output. Therefore, to get a set of different perturbations along a given single trajectory, a set of 
different starting random numbers should be used. Note, however, that if any other parameters are 
changed (different spacing along the trajectory, different starting position, etc.) then the same start- 
ing random number will produce a different set of random perturbations. 

6. 4-D Scratch Unit, IOTEM1: In order to save array space, the 4-D profiles required to 
interpolate to the 5° x 5° grid locations are read from the 4-D files to this scratch file rather than 
being put into arrays. The unit number for this scratch file can be any available unit. Normally, 
the file is a temporary (scratch) file, and, if so, does not have to be assigned before execution of 

the program. 

7. NMC Grid Point Scratch Unit, IOTEM2: Also in order to save computer storage, the 
NMC grid point array read in from the SCIDAT file is stored in a temporary scratch file and does 
not have to be assigned before execution of the program. 

8. 4-D Diagnostic Print Output Option, IPRT: Use a value of zero to get diagnostic 
information printed out at the initial generation of a 4-D grid. Use IPRT = 1 to suppress this 
diagnostic output. 

9. Maximum Allowable Test Violations, NLIMIT: In the generation of 4-D data grids, 
GRAM-90 tests for perfect gas law and hydrostatic law violations in the original 4-D data base. 
NLIMIT sets the upper limit on the allowed number of test violations which can occur without 
automatic termination of the program. Normally, NLIMIT = 12 is used, but a value up to 
NLIMIT = 25 is allowed. Any input value of NLIMIT less than 12 is reset to 12; any value 

greater than 25 is reset to 25. 

Input Line Number 3. This card is read by the SETUP subroutine and contains starting 
values for the random perturbation parameters at the initial position and the scale parameter for the 
random perturbation magnitudes. If random perturbations are not to be computed (Random Option 
- 2), this line should not be put in. All values of this free-field format line are real. For a normal 
design application, the values of the initial perturbations on this line should all be zero, unless the 
run is to be a continuation of a previously run trajectory or profile segment, in which case the 
output random parameters of the last output position are input, and the last output position becomes 
the initial position of the new run. 

1-6. Initial Pressure, Density, and Temperature Perturbations: Parameters RP1L, 

RP1S, RD1L, RDIS, RT1L, and RT1S are, respectively, the initial values of random relative 
pressure (p'/p), density (p'/p), and temperature (T'/T), in percent, for the large-scale (L) and 
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small-scale (S) components. These are starting values for the initial position. Use zeros for standard 
design applications. 

7-10. Initial Wind Component Perturbations: Parameters RU1L, RU1S, RVIL, and 
RV1S are initial values of the random eastward (U) and northward (V) random wind components 
in m/s for the large-scale (L) and small-scale (S) components. Use zeros for standard design 
applications. 

1 1 . Random Perturbation Scale (RPSCALE): This is a parameter to allow random per- 
turbations which are more disturbed (RPSCALE > 1) or more quiescent (RPSCALE < 1) than 
normal (RPSCALE = 1). Values can range from zero (no perturbations) to two (twice normal 
perturbation magnitudes). 

Trajectory Input. The free-field trajectory position input and backup record are put in 
only if a trajectory is to be evaluated, rather than a linear profile, generated automatically in the 
program from information on the first input line. There is no limit to the number of trajectory posi- 
tion records which can be put in. The program continues evaluating the atmospheric parameters and 
looping back to read a new trajectory position until a position below the surface is reached, or 
until the trajectory backup record is reached. Each free-field trajectory record has the time (integer 
seconds), the height (kilometers), the latitude (degrees, southern latitude negative), and the west 
longitude (degrees, 0° to 360° or east longitudes negative). Any east longitudes read in as negative 
values are converted to the 0° to 360° system before being used by the program. The trajectory 
backup record has the same free-field form as a regular trajectory record, except that any negative 
value for height is used. The negative height terminates the loop which evaluates atmospheric 
parameters. The trajectory data can be input from any unit (with trajectory option = unit number). 
The trajectory option is item 18 on line number 1. 

4.5 Output of the Program 

The first few lines of print output are primarily a listing of the input parameters. Following 
a heading which describes each output value for the trajectory or profile evaluations, the position, 
time monthly mean, and total pressure, density, temperature, and winds are listed for each posi- 
tion. The thermal wind shear for the monthly mean winds, the percent deviation from the standard 
atmosphere, the mean vertical wind, and the perturbation data are also given for each position. The 
perturbation data consists of the stationary perturbations, the quasi-biennial values at the position 
and time, the quasi-biennial magnitudes, the random perturbation values, and the random perturba- 
tion standard deviations. Optional nonprint (e.g., disk) output for values at each position is also 
available to be used for input to plotter programs, or for other purposes. 

Heading Information. Primarily, the heading information contains a listing of the input 
data values. However, there are some changes from the values input. If an east longitude is put in 
as a negative value, — 180° < lat. < 0°, then it is converted to a west longitude in the 0° to 360° 
range before the heading is listed. The program evaluates the initial random pressure, density, 
temperature and wind standard deviations and the initial density-velocity correlation from data in 
the SCIDAT9 data file, and lists the computed values on the heading. The Julian date is computed 
by the program from the input date and is also listed with the heading information. The Julian date 
is required by the Jacchia and QBO sections of the program. 
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Position and Time Output. Positions and times as generated by the automatic linear 
profile features, or as input by the trajectory input data, are listed on the output. The time is given 
in seconds. Within the program, the input time in hours, minutes, and seconds are updated in that 
form also. However, only a continuously increasing time in seconds is printed out. If time in 
hours, minutes, and seconds were desired, these variables could easily be printed^out by adding 
them to the output list. All output west longitudes are converted to the 0° to 360 range before 
being printed out. If a latitude greater than 90° in absolute magnitude is generated (or input) then a 
transformation is made. 
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long. = long. + 180° . 

(4.7) 


Monthly Mean Data. The monthly mean values of pressure, density, and temperatures 
consist of either: (I) values from the 4-D data file if the height is below 25 km, (2) the sum of 
middle atmosphere zonal-mean plus stationary perturbation values if the height is between 30 and 
90 km, (3) an interpolation between 4-D at 25 km and zonal-mean plus stationary perturbations at 
30 km’ if the height is between 25 and 30 km, (4) Jacchia (MET) model values if the height is 
above 120 km, or (5) faired values between middle atmosphere and MET model values d the 
height is between 90 and 120 km. 

The percent deviations from the U.S. 1976 Standard Atmosphere are evaluated by using 
standard atmosphere values computed by the subroutine STD ATM. The percent deviations are eval- 
uated by the relations 100 (T-T S )!T S , 100(p- p,)/p„ and \00(p-p s )p s , where the subscript 5 refers 
to the standard atmosphere values. This subroutine accurately reproduces the tabulated U.S. 1976 
Standard Atmosphere values to within an accuracy of better than 0.2 percent above 90 km. The 
subroutine reproduces the tabular values even more accurately in the height region below 90 km, 
where the molecular weight is constant. Since the U.S. 1976 Standard Atmosphere is not defined 
above 1,000 km, the percent deviations printed out for heights above 1,000 km are zero. Since the 
Jacchia (MET) model is sensitive to solar activity conditions, large deviations from U.S. Standard 
Atmosphere values can be produced in this height range for certain ranges of F10.7 and a p values. 

The thermal wind shear values are values of du/dz and dv/dz for the monthly mean wind 
(see section 2). The wind values are computed from the geostrophic wind equation or from the 
new wind data base. The wind shear components are computed by either the thermal wind 
equations, which are determined by the horizontal gradients of the monthly mean temperature, or 
computed as vertical gradients from the new wind data base. The mean vertical wind is evaluated, 
as described in section 2, by combinations of horizontal and vertical temperature gradients and the 
monthly mean winds. 
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The Total (Mean Plus Perturbation) Data. The parameter values listed under the head- 
ing of “Mean Plus Perturbations” are the monthly mean values, as defined above, plus the random 
perturbations, plus (if the height is between 10 and 90 km) the quasi-biennial perturbations. These 
mean-plus-perturbation values represent values which would be typical “instantaneous” values of the 
pressure, density, temperature, and winds. The percent deviations from the U.S. Standard 
Atmosphere are computed in the same way as for the percent deviations of the monthly mean 
values from the standard atmosphere. 

Perturbation Values. The data under the “Perturbation Values" heading are the various 
perturbation values, magnitudes, and amplitudes. The stationary perturbations (denoted SP on the 
printout) are defined only if the height is between 25 and 90 km. The monthly mean, y,„, of 
parameter v should be the zonal-mean value Z y , evaluated from the SCIDAT9 data file, modified 
by the given stational perturbation value s y , in percent, by the relation 


y,„ = Z y ( 1 +yi00) . 


(4.8) 


The data labeled “QBO” are the values of the QBO oscillation at the output time and position. The 
data labeled “MAG” give the magnitude of the QBO oscillations at the output position and time. 
The QBO perturbation values should always be less than or equal to the magnitude values in 
absolute value. The data labeled “RANL,” “RANS,” and “RANT” are the large-scale, small-scale, 
and total-random perturbations evaluated at the output time and place. The data labeled “SIGL,” 
“SIGS," and “SIGT” are the standard deviations of the large-scale, small-scale, and total-random 
components at the output time and positions. According to the Gaussian distribution, on which the 
random perturbations are based, the perturbation values should be within the range ±o 68 percent 
of the time and outside the range ±a 32 percent of the time. Similarly, the perturbation values 
should be within the range ± 2a 95 percent of the time and outside the range ± 2a 5 percent of 
the time. The evaluation of the QBO and random perturbation output can be suppressed by the 
QBO and random options, if desired. SIGW gives the magnitude of the random vertical wind 
perturbations (m/s). These data are from the turbulence model of Justus et al. (1990). 

Nonprint Output. The nonprint output is available, as an option, controlled by the input 
value of the output option parameter. If nonprint output is desired, it comes out in the form of 
records with format F5.1, 7F7.2, 4F6.2, containing the following information: (I) the height in 
kilometers; (2) the latitude in degrees; (3) the west longitude in 0° to 360°; (4) the percentage devi- 
ation of the mean monthly value of density from the 1976 U.S. Standard Atmosphere; (5) the 
monthly mean temperature; (6-8) the eastward, northward, and vertical components of the monthly 
mean (geostrophic) wind; and (9-12) the magnitudes of the large-scale random perturbations in 
pressure, density, temperature (percent), and eastward and northward wind (m/s). If different output 
parameters or formats are desired, the statement WRITE (IOPP, 951)..., immediately after label 
920, and/or the format 951, in subroutine SCIMOD, can easily be modified to accommodate these 
requirements. 


48 



4.6 Program Diagnostics 

There are several possible reasons which can cause the printing ot diagnostic messages and 
termination of the run during the setup phase. If during the setup procedure, the NMC grid point 
number data table does not contain the required 1 ,977 values, a message N RECORDS 
WRITTEN BY SETNMC IN SCRATCH FILE M” is printed, and execution is terminated. A 
premature end-of-file encountered while attempting to read the NMC data from the SCIDAT9 file 
results in a diagnostic message “PREMATURE END-OF-FILE FOUND ON UNIT IUS CALLED 
FROM SUBROUTINE GETNMC” where IUS is the SCIDAT9 unit number. A similar message is 
produced by subroutine RTRAN if a premature end-of-file occurs elsewhere in the SCIDAT9 file. 

If, during the reading of the SCIDAT9 data file, any record is read which does not have the 
expected code character or characters, then the message results: “ERROR IN SETUP INPUT 
followed by a listing of the latest data values read in. This message is also produced it the random 
option and the quasi-biennial option do not have a value of either one or two. Any condition which 
results in this error message terminates the execution. 

If an initial random number outside the range 1 to 30,000 is input, the message FIRST 
RANDOM NUMBER OUT OF RANGE” is printed and execution is terminated. 

There are also several conditions which could result in diagnostic messages in the 4-D sec- 
tion: If, during the reading of the 4-D data file, the incorrect month or record number is encounte- 
red then a message “***** UNIT NO. IU4 IN ERROR. IRN = XXXX DTEMP(213) - XXXX 
MONTH = XX DTEMP(212) - XXXX IPT(I, J) = XXXX” is printed, where: 

IU4 = Unit on which 4-D file is being read 


IRN = 4-D file record number 

DTEMP(213) = Month value in last record read from 4-D file 


MONTH = Month requested to be read 

DTEMP(212) = Point number in last record read from 4-D file 

IPT(I, J) = Point number required for profile J to be interpolated to Ith requested profile. 


If this error condition occurs, execution is terminated. 


If, during steps along a profile or trajectory, the next position is outside the current 4-D 
grid area! a new 4-D grid is generated, suitable for the new location. If the new ^ 
successfully generated, a message “UNABLE TO GENERATE 4-D GRID. TOO MANY RETRIES 
IN INTER4” is generated and program execution terminates. 


If a large number of perfect gas law violations are encountered or a hydrostatic violation 
occurs at too low an altitude, the program may terminate execution with a message “UNABLE TO 
GENERATE 4-D GRID. TOO MANY TEST VIOLATIONS,” or “UNABLE TO GENERATE 
4-D GRID. IHV = XX,” where IHV is the height where the lowest hydrostatic violation occurred. 
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In some cases, this condition may be removed by increasing the value input tor NLIMIT (item 9, 
input line 2). 

As optional diagnostic output, various information can be requested (by parameter IPRT: 
item 8, input line 2). These output are generated in the 4-D section of the program, during the 
generation of the 4-D grid, by subroutines GEN4D, GRID4D, and SCIMOD, This information 
includes data on which 4-D data produced gas law violations or hydrostatic law violations, which 
4-D values had to be replaced by new values as a result of the gas law and hydrostatic tests, along 
with the new values generated. The complete array of 4-D pressure, density, and temperature data 
(means and standard deviations) is also printed out, both those read in directly from the 4-D file, 
and those produced after the gas law and hydrostatic tests. 

5. PROGRAMMER’S MANUAL 
5.1 Description of Subroutines 

The following is a brief description of each of the GRAM-90 program subroutines, in 
alphabetical order: 

ADJUST: Adjusts the 4-D profiles of pressure, density, and temperature variance (read from 

the 4-D file) to satisfy the Buell constraints imposed by the perfect gas law. 

CHECK: A consistency check routine for the 4-D grid data produced by GEN4D. CHECK is 

called for each height to be evaluated, and tests for reasonable values of scale height 
immediately above and below that height. It also tests for reasonable horizontal 
pressure gradients. Failure of either test produces the diagnostic asterisk between the 
output values of wind components. 

CORLAT: Evaluates the horizontal and vertical scales for large- and small-scale density, 

temperature, and wind components; computes the autocorrelations and cross correla- 
tions for the two scale perturbation model; and evaluates new perturbation values 
having appropriate correlations with the perturbations at the previous position. 

CORREL: Correlation function used in the random perturbation model. 

DIAGEQ: A matrix diagonalizing procedure used by the ADJUST subroutine. 

FAIR; Fairs between the middle atmosphere and Jacchia (MET) values in the 90- to 120-km 

height range. 

FAIR5: Fairs between the region below and above 500 km in the MET subroutine. 

GAUSS: Gaussian quadrature routine, used in the MET model. 

GEN4D: Generates the 4-D polar ( | latitude) > 75°) or nonpolar (16, 5°x5° points) grid of 

pressure, density, temperature, and variance profiles. 
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GETNMC: 

GRAM: 

GRID4D: 

GROUP: 

GTERP: 

INTER W: 

INTERZ: 

INTER2: 

INTER4: 

INTLL: 

INTRP4: 

INTRUV: 


Reads the NMC grid point values from the SCIDAT9 data file and loads them onto 
a scratch file. This subroutine is essentially unchanged from the subroutine of the 
same name in the original 4-D program. 

The main segment of GRAM-90. The main segment serves as a driver for the 
SETUP and SCIMOD subroutines. 

After the array of 4-D grid latitudes-longitudes has been evaluated, this subroutine 
looks up the data from the 4-D data file and interpolates to determine profiles of 
pressure, density, temperature, and variance at the 4-D grid locations. Protiles to be 
interpolated to 4-D grid locations are loaded onto a scratch file from the 4-D file 
before the interpolation is done. 

A subroutine, called by CHECK, which groups the 16 4-D pressure data at the given 
height into one or more groups which have consistent and reasonable horizontal 
pressure gradients within each group. If the subsequent geostrophic wind calculations 
in WIND use horizontal pressure gradients evaluated from differences across incon- 
sistent groups of 4-D data, the diagnostic asterisk is printed between the output 
values of wind components. 

Uses linear latitude interpolation, linear temperature, and linear logarithm of density 
interpolation on height to evaluate middle-atmosphere zonal-mean data to a given 
latitude and height. 

Two-variable, linear interpolation between known value U1 and VI at height Z1 and 
U2 and V2 at height Z2 to determine U and V at height Z, where Z is between Z1 
and Z2. 

Three-variable interplation, linear on temperature, and gas constant ( R = p/p71. and 
linear on the logarithm of density, with pressure computed from perfect gas law and 
interpolated temperature, density, and gas constant. 

Three-variable interpolation, linear on all three variables. 

Interpolates between the pressure, density, and temperature profiles at the 4-D grid 
locations. This subroutine calls subroutine INTLL to do the latitude-longitude 
interpolation. 

One-variable interpolation between values in an array of latitude and longitude loca- 
tions by equation (2.4). 

The subroutine for the latitude-longitude interpolation of values from the 4-D data 
file into the 4-D grid array. This is a modification of the INTERP subroutine of the 
original 4-D program. 

Evaluates the standard deviations of the random wind components at given height 
and latitude by calling INTERW subroutine. 
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INTRW: 

INTR25: 

JAC: 

JACCH: 

J70: 

MOLWT: 

PDTUV: 

PERTRB: 

PHASE: 

PPND: 

QBOGEN: 

RANDOM: 

RIG: 

RTERP: 

RTRAN: 


Finds standard deviation for random vertical wind components. 

Subroutine to find large-scale fractional variances (f L in equations (4.2) and (4.3)) at 
a given height and latitude. 

Calculates the molecular weight, density, and temperature for the Jacchia (MET) 
model. 

Main subroutine driver for the Jacchia (MET) section, serves as a driver for J70 and 
other Jacchia section subroutines. 

Subroutine called by JACCH driver subroutine for the Jacchia (MET) model. 
Calculates molecular weight for the Jacchia (MET) model. 

Interpolates the stationary perturbations on latitude and longitude at a given height. 
This subroutine is similar to INTLL. 

Evaluates the pressure, density, temperature, and wind component random perturba- 
tions by the correlated random perturbation model discussed in section 2. 

A linear height-latitude interpolation routine for the quasi-biennial phase. The 
interpolation properly accounts for the phase discontinuity between 0 and 870 days 
(the quasi-biennial period). 

Function to produce normal (Gaussian) distributed random variable, for use in the 
random perturbation model (q t , in equation (2.19)). 

Computes the QBO perturbation values and their amplitudes and phases. The ampli- 
tudes and phases of the QBO pressure, density, temperature, and wind perturbations 
are interpolated from the amplitude and phase data from the SCIDAT9 data file, by 
calling the INTERZ and INTERW subroutines. 

Produces a random number selected from a uniform distribution between zero and 
one. This is required as input to the subroutine PPND. 

Computes the acceleration of gravity and the radius from the center of the Earth for 
a position at a given latitude and height. 

Computes the standard deviations of the random pressure, density, and temperature 
perturbations by calling subroutine INTERZ. 

This subroutine contains several read sections with multiple entry points coming from 
subroutine SETUP. The read statements are for reading the SCIDAT9 data file. 
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SCIMOD The heart of the GRAM-90 program. This subroutine branches on height to evaluate 
the atmospheric parameters by the Jacchia (MET), the middle atmosphere, or the 4-D 
methods The QBO and random perturbations are also evaluated and the output is 
printed (and optionally also written to a file) by the SCIMOD subroutine. 

SELEC4: Selects the 4-D data needed for interpolation. This subroutine is a modification of 

the INPUT subroutine of the original 4-D program. 


SETUP: 


SORT4: 


This subroutine reads in the NMC grid points with the GETNMC subroutine and 
reads and loads the data from the required month in the SCIDAT9 data file into 

arrays. 

Computes the seasonal-latitudinal variation of density in the lower thermosphere por- 
tion of the Jacchia (MET) model. 

Sorts the 4-D locations for sequential reading from the 4-D data file. This subroutine 
is a modification of the SORT subroutine from the original 4-D program. 


STD ATM: Evaluates the 1976 U.S. Standard Atmosphere values of pressure, density, and 

temperture at any given height up to 1 ,000 km. 


TINF: This subroutine computes the exospheric temperature for the Jacchia (MET) model. 


This subroutine calculates the variables necessary for input into the subroutine TINF 
in the Jacchia (MET) model. 


WIND: This subroutine evaluates the geostrophic winds from input values of horizontal 

pressure gradient if the height is less than 25 km or more than 90 km. If the latitude 
is below the minimum geostrophic latitude, it evaluates geostrophic wind at 
minimum geostrophic north latitude and at minimum geostrophic south latitude and 
then interpolates in between. If the height is between 25 and 90 km, the new middle 
atmosphere wind data base is used. Between 20 and 25 km and between 90 and 95 
km, a smooth fairing between geostrophic and middle atmosphere wind is used. 


5.2 The Primary Section 

This section consists of the main program segment GRAM, the SCIMOD subroutine the 
subroutines for evaluating the middle atmosphere values, the stationary perturbations, the QBO and 
random perturbations, and general interpolation subroutines. With the exception of GRAM and 
SCIMOD, the parts of this section were adequately described in the previous section. 

Many of the subroutines transfer their input and output via COMMON statements. This 
procedure saves much in core storage space. 

Main Segment GRAM. This program serves as a driver for the SETUP and SCIMOD 
subroutines. It reads one line of input, the first input line, as described in section 4. The trajectory 
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input records (if used) are also read by GRAM, after control has returned from SETUP, which 
reads the second and third initial data input lines. Section 4.4 and appendix B give a further 
description of the input required. 

Subroutine SCIMOD. This is the primary subroutine of the GRAM program. It serves as 
a driver for all of the various sections of the atmospheric evaluation. The SCIMOD subroutine 
prints and (optionally) generates output on a nonprint output file, the output described in section 4 
and appendix B. It also transfers data values to other subroutines via the above-mentioned 
COMMON lists. The SCIMOD subroutine updates the profile or trajectory positions by setting the 
current position equal to the previous position before exit. The previous position information then 
stays in the COMMON lists until the next call to SCIMOD. The previous random perturbations are 
handled in similar fashion. 


5.3 The Setup Section 

The function of the setup section of the program is to load the initial data and the data from 
the SCIDAT9 file. 

The SETUP subroutine reads the second and third lines of input, as described in section 4. 
COMMON lists transfer the arrays, loaded with the appropriate data from the SCIDAT9 data file, 
to the other subroutines. COMMON lists are also used to transfer data to SETUP from the read 
subroutine RTRAN, which has multiple entry points, for various different types of data from the 
SCIDAT9 data file. 


5.4 The Jacchia Section 

The subroutine JACCH calculates the pressure, density, and temperature at a point in space 
for heights above 90 km for a particular time. It is essentially unchanged from the MET model 
(Hickey, 1988a,b). 


5.5 The 4-D Section 

GRID4D and subroutines SORT4, INTRP4, and SELEC4 are basically the MAIN 
PROGRAM, SORT, INTERP, and INPUT as documented in the 4-D user’s reference manual and 
subsequent updates (Fowler and Willard, 1972; Spiegler and Fowler, 1972). Some changes have 
been made, as documented in earlier versions of the GRAM program (e.g., Justus et al., 1980). In 
GRAM-90 a new option has been added which allows output of extensive diagnostic data from the 
4-D section during the generation process for producing the 4-D grid values. 

5.6 GRAM-90 Input and Output 

Input required by GRAM-90 for setting up the desired profile or trajectory, and for selecting 
the various program options, was described in section 4. Appendix A gives more complete details 
on the middle atmosphere data file, SCIDAT9 GRAM-90 output was also discussed in section 4. 
Additional details on program input and output are given in appendix B. 
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APPENDIX A 


DESCRIPTION OF THE REVISED MIDDLE ATMOSPHERE DATA FILE 
“SCIDAT9” FOR THE GRAM-90 PROGRAM 


The SCIDAT9 file contains the following data, identified by code characters at the 
beginning of each record. For additional discussion of the SCIDAT9 input data, see section 4.3. 


Code 

Data 

Description 

N 

NMC grid data 

Same as NMC grid required by original 4-D 
program. Data consists of sequential point 
number followed by the two corresponding 
NMC grid indices. There are five points 
per record in format (A2, 1517). 

ZP 

ZD 

ZT 

ZU 

Zonal mean pressure (N/m 2 ) 
Zonal mean density (kg/m 3 ) 
Zonal mean temperature (K) 
Zonal mean zonal wind (m/s) 

Month, height, values at latitudes 
-90, -80, ..., 10, 20, ... 90, 
exponent. Format is A2, 14, 15, 1916, 
14. 

SP 

Stationary perturbations in monthly 
mean pressure (per mil) 

Month, height, latitude, A p/p at longi- 
tude 180, 160W, I40W, ..., 140E, 160E, 
in format A2, 2115. 

SD 

Stationary perturbations in monthly 
mean density (per mil) 

Same as for code SP data, for Ap/p. 

ST 

Stationary perturbations in monthly 
mean temperature (per mil) 

Same as code SP data, for A777\ 

su 

Stationary perturbations in monthly 
mean zonal (eastward) wind component 
(0.1 m/s) 

Same as for code SP data, for Au. 

sv 

Stationary perturbations in monthly 
mean meridional (northward) wind 
component (0.1 m/s) 

Same as for code SP data, for Av, 

RP 

Random perturbation magnitude for 
pressure (percent) 

Month, height, values of a p /p at latitudes 
-90°, -80°, ..., 80°, 90°, exponent, 
in format A2, 14, 15, 1916, 14. 

RD 

Random perturbation magnitude for 
density (percent) 

Same as for code RP data, for a>/p. 
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Code 


Data 


Description 


RT 

RU 

RV 

P 

PW 

CS 

CL 

QP 

QD 

QT 

QU 

QV 


Random perturbation magnitudes for 
temperature (percent) 

Random perturbation magnitudes for 
zonal wind (m/s) 

Random perturbation magnitudes for 
meridional wind (m/s) 

Fractional variance in large scale 
thermodynamic variables 


Fractional variance in large scale wind 


Small-scale density-velocity correlations 


Large-scale density-velocity correlations 


QBO pressure parameters-amplitude 
(per mil) and phase (days after January 
0, 1966 when first maximum occurs) 

QBO density parameters (as in QP) 

QBO temperature parameters 

QBO eastward wind parameters-amplitude 
(0. 1 m/s) and phase (days after 
January 0, 1966) 

QBO northward wind parameters (as in 
QU) 


Same as for code RP data, for oy/7’. 


Same as for code RP data, for a„. 


Same as for code RP data, for ct v . 

13 (denotes annual average) height, 
fractional variance in large scale per mil 
for pressure, density, and temperature, each 
at latitude 10°, 30°, 50°, 79°, 90°, 
in format A2, 1717. 

13 (denotes annual average) height, 
fractional variance in u at 10°, 30°, 50°, 70°, 
90°, in format A2, 1217. 

13 (denoting annual average), height, 

<pw>s at 10°, 3(F, 50°, 70°, 90° latitude, 
same for <pv> s , in format A2, 1217. 

13 (denoting annual average), height, 

<P«>/. at 10°, 30°, 50°, 70°, 90° latitude, 
same for <pv>/, in format A2, 1217. 

Height, amplitude, and phase at 10° 
latitude, amplitude, and phase at 30° ..., 
amplitude and phase at 90°, in format 
A2, 1117. 

Same as for code QP data. 

Same as for code QP data. 

Same as for code QP data. 


Same as for code QP data. 
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The SCIDAT9 file contains the NMC grid data, the middle atmosphere zonal mean data and 
stationary perturbation data, the random perturbation data, the fractional large-scale variances and 
t^e densTty^ -velocity correlations, and the QBO data. Each record of the NMC grid data fde con- 
tains the code (N) and .v-v coordinates for five points. The total number ol N gri poin . 

W7 The NMC grid data file contains a total of 396 records, with the last record containing points 
1 976 and 1 977 and zeros for the remaining values. The zonal mean data contains , recor . , 
the!nattonary perturbation data contains 15,300 records, the code random perturba, ton data con a, ns 
1 596 records the code P large-scale fractional variances contain 25 records, the code PW g 
scale* fractional w!nd variance! contain 25 records, code CS and CL density-veloctty correlate 
data contain 25 records each, and the QBO data contain 80 records. 
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Line 


APPENDIX B 


SAMPLE INPUT AND OUTPUT FOR THE GRAM-90 PROGRAM 


Input to GRAM-90 is as follows: 

(All input data lines are in free-field format. R designates real; I denotes integer). 


INITIAL HEIGHT (R) 


Height of starting position, km. 


INITIAL LATITUDE (R) 


INITIAL WEST LONGITUDE (R) 


10.7 cm FLUX (R) 


MEAN 10.7 cm FLUX (R) 


GEOMAGNETIC INDEX (R) 


DATE: MONTH, DAY, YEAR (I) 


GREENWICH TIME: HOUR, 

MINUTE, SECOND (I) 


LAT INCREMENT (R) 


Latitude of starting position (degrees, southern latitudes 
negative). 

West longitude of starting position (degrees, 0° to 360°, 
or east longitudes negative). 

Solar 10.7 cm radio noise flux (10 22 watts/m 2 ) at time 
of calculations. Use zero if height does not go over 90 km. 
Use 230 for design applications or consult MSFC for monthly 
predictions. 

81 day mean solar cm flux. Use zero if height does not go 
over 90 km. Use 230 for design applications or consult 
MSFC for monthly predictions. 

Geomagnetic index a p . Use zero if height does not go 
over 90 km. Use 20.3 for design steady-state conditions, or 
400 for maximum conditions, or consult MSFC. 

Date for starting time of calculations (month, date, two 
digit year) 

Time for starting position (hours, minutes, seconds). 

Use time corresponding to local time = 0900 for design 
steady-state, or local time = 1400 for maximum conditions. 

Latitude displacement (degrees) between successive positions 
(new lat. = old lat. + lat. increment). Use zero if trajectory 
positions are to be read in. 


WEST LON INCREMENT (R) West longitude displacement (degrees) between successive 

positions (new long. — old long. + long, increment). Use zero 
if trajectory positions are to be read in. 

HEIGHT INCREMENT (R) Height decrease (km) between successive positions 

(new height = old height - height increment). 

Normal profiles are generated downward. If an upward 
generated profile is desired, set height increment negative. 
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Line 2 ► i Line 1 (cont.) 


MAXIMUM NUMBER OF 
POSITIONS (I) 

TIME INCREMENT (I) 

TRAJECTORY OPTION (I) 

OUTPUT OPTION (I) 

MIN. GEOSTROPH. LAT (R) 


SCIDAT9 INPUT UNIT (I) 

4-D INPUT UNIT (I) 

RANDOM OPTION (I) 

QBO OPTION (I) 

FIRST RANDOM NUMBER (I) 

4-D SCRATCH UNIT (I) 

NMC GRID POINTS SCRATCH 
UNIT (I) 

4-D DIAGNOSTIC OUTPUT 
OPTION (I) 

MAX ALLOWABLE TEST 
VIOLATIONS (I) 


Number of positions to be computed, not including initial 
position. Use zero if trajectory positions are to be read in. 

Time displacement (seconds) between successive positions for 
automatically generated profiles (new time = old time + 
time increment). 

0 for linear profile generated automatically, internal to the 
program, or value equal to unit number for a trajectory with 
each position to read in. 

0 for no nonprint output of atmospheric parameter values, 
or value equal to unit number to get nonprint output. 

Lowest latitude (magnitude) for which geostrophic winds are 
to be used in 4-D (0-25 km) and Jacchia (above 90-km height 
segments). Otherwise, interpolation is used to fill in winds. 

In middle heights (25-90 km), the new middle atmosphere data 
base is used at all latitudes. 

Unit number for file containing SCIDAT9 data. Use any 
available unit number. 

Unit number for 4-D input data file. Use any available 
unit number. 

Use 1 to compute random perturbation output, 2 means 
do not compute random perturbation output. 

Use 1 to compute QBO output, 2 means do not compute 
QBO output. 

Initial number for random generator used to compute random 
perturbations (can be any positive integer 1-30,000). Use I for 
standard design applications. 

Unit number for scratch file for 4-D grid profiles required 
in computations. Use any available unit number. This normally 
is a temporary scratch file. 

Unit number for scratch file to store NMC grid point data. 

Use any available unit number. This normally is a temporary 
scratch file. 

Use 0 for optional 4-D diagnostic output, and 1 
for no diagnostic output. 

Maximum number of test violations allowed for 4-D data, 
normally 12, can be up to 25. 
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ii INITIAL PL, PS, DL, DS, Initial values of large-scale and small-scale random relative 

TL, TS (R) pressure, density, and temperature perturbations, percent. Use 

zeros for standard design applications. 

g INITIAL UL, US, VL, VS (R) Initial values of large-scale and small-scale random wind 

j components, m/s. Use zeros for standard design applications. 

RPSCALE Random perturbation scale: 1 = normal, minimum 0, 

jf maximum 2. 


^Include line 3 only if random option = 
TRAJECTORY INPUT (II, 3R) 



Use only if linear profile is not to be generated 
automatically. Each record has time (seconds as integer), 
height (km), latitude (degrees), and west longitude (degrees). 


TRAJECTORY BACKUP RECORD Only if trajectory input is used. Same form as a trajectory 
(II, 3R) position but with any negative height value. 


The trajectory input records are optional, in free-field format. If included, use as many 
records as necessary. 

Input for the following sample output listing is as follows: 

LINEl: 142.0, 28.45, 80.53, 230.0, 230.0, 20.3, 1, 1 , 90, 0, 0, 0, .0, .0, 2., 72, 0, 0, 99, 18 

L1NE2: 3, 4, 1, 1, 1, 12, 13, 1, 12 

LINE3: 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1. 

A SUMMARY OF THE ORGANIZATION OF AN INPUT DATA FILE IS AS FOLLOWS: 

Initial Data 

Line I, as described at the beginning of this appendix. 

Line 2, as described at the beginning of this appendix. 

Line 3, optional, included only if random option = 1. 

Trajectory Position Data and Backup Card 

Optional. Include if trajectory, rather than linear profile generated by the program is to be 
evaluated. 
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More Data of the Same Kind (Starting With Initial Data, Line I) 

If additional trajectories or profiles are to be evaluated, the data may be input one set 
immediately after the other. The program is actually more efficient for such multiple runs if the 
month remains the same. This is because as long as the month remains the same, the SCIDAT9 
data file read can be avoided for each subsequent data set. 

OUTPUT TO GRAM: 

INITIAL POSITION Listing of input values of initial position. 

PROGRAM OPTIONS Listing of input values of program options. 

HEIGHT, LAT., LONG., TIME Position and time where atmospheric parameters are evaluated. 

UNPERTURBED PRESSURE Computed from Jacchia (MET), 4-D, or middle atmosphere 
DENSITY, TEMPERATURE, (zonal mean plus stationary perturbations), depending on height. 
AND WIND (monthly mean 
values) 

TOTAL PRESSURE, Monthly means plus random perturbations and QBO 

DENSITY, TEMPERATURE, perturbations. 

AND WIND 

THERMAL WIND SHEAR From thermal wind equations using finite differences of 

Jacchia (MET), 4-D, or middle atmosphere data, depending 
on height. 

MEAN VERTICAL WIND From mean isentropic surface slopes. 

PERTURBATION VALUES Stationary perturbations, QBO perturbations, and amplitudes, 

and random perturbations and magnitudes for the small-scale 
(S), large-scale (L), and total (T) perturbations. Perturbations 
are those which are added to monthly means to produce total 
results output. 

Following is a listing of sample output from the GRAM-90 program. Initial lines of output 
are merely listings of the input data for each reference. For a listing of the input data for this 
sample output, see above. 
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ERRATA 


THE NASA MSFC GLOBAL REFERENCE ATMOSPHERIC 
MODEL— 1990 VERSION (GRAM-90) 

PART II FICHE ERRATA 

l sers ol the GRAM-90 part II program and data listings, from the enclosed Fiche. will 
need to apply the following errata corrections in order to interpret a correct listing. 

FicheOi-Row 3-Column 6-Page 29 

In SUBROUTINE INTRUV- 

Imertline INUV ?A to read: 


I=l+INT(Hy5 INUV 6 

IF (H.GE. 125) I=25+(INT(H)-i20)/29 INUV 7 

IF (I.GT. 29) 1=29 INUV 7A 

C UPPER HEIGHT INDEX INUV 8 

IP=I+1 INUV 9 


FicheOl-Row 4-Column 8-Page 45 
In SUBROUTINE RTERP- 

Inserthne RTER 9A and modify fetes RTER 11 and RTER 25 to read: 


IF(H GE 120.)I=25+INT((H-12O.y2O.) RTER 8 

IF(I.LT.1)I = 1 RTER 9 

IF(I.GT.29JI = 29 RTER 9A 

IP = 1+1 RTER. 10 

IF(IP .GT. 29)IP=29 RTER 11 

C LOWER LATITUDE INDEX RTER 12 

J=INT((PHI+1S0.V10.) RTER 13 


72=5 *ip-5. RTER 23 

GO TO 40 RTER 24 

30 Z2=l2O.+20*(lP-25) RTER 25 

40 PHI1=-100.+1O.*J RTER 26 

PHI2=-100+10.*JP RTER 27 
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Me 81 - Row 4 - Column 12- Page 49 


In SUBROUTINE SCIMOD- 

Modify lines 3CIM 141, SCIM 158 and SCIM 159toiead. 

THETE=THET-5. SCIM 140 

C.....JACCHIA VALUES AT LOVER HEIGHT. CURRENT LON-LAT+5 DEGREES SCIM 14 1 
C UT, FOR DP/DY AND DT/DY SCIM 142 


PHIN=PHIR+5./TAC SCIM 157 

THETE=THET-5. SCIM 158 

..JACCHIA VALUES AT UPPER HEIGHT, CURRENT LON-LAT+5 DEGREES SCIM 159 

LAT, FOR DP/DY AND DT/DY SQM 160 
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